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1 Introduction 

A possible approach to explaining the observational evidence for cosmic acceleration is that 
general relativity is modified in regions of low density, low acceleration or, simply, on large 
scales. A number of models have been proposed [1], and the idea that there may be a gravita¬ 
tional solution to the dark energy problem has led to a renewed scrutiny of the fundamental 
properties of gravity. The hope is that, in the very least, a better understanding of general 
relativity will emerge from current and future cosmological observations. 

A common feature in a number of modified theories of gravity is that there is some form 
of screening [2]. Precision tests [3, 4] indicate that general relativity is at least an excellent 
effective theory for describing gravitational physics locally or in regions of moderately high 
density and curvature. This means that, in any modification of the theory, the nature 
of gravity must depend on its environment and that deviations from general relativity are 
absent (or screened) in regions in which it has been well tested. Gurrent proposals involve 
the chameleon [5, 6], symmetron [7-9] or Vainshtein [10] mechanism, all of which lead to 
particular features in the gravitational potential. 

A key characteristic of the various screening mechanisms [5, 7, 10] currently being 
considered is that they are fundamentally non-linear. While they might have an effect in 
regimes where density perturbations are still in the linear regime, their full effect comes into 
play in the non-linear regime of gravitational collapse. This means that, if we are to fully 
understand the effect of screening in universes which closely resemble our own, we must be 
able to model non-linear structure formation not only accurately but also efficiently. There 
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has been substantial progress in numerically modelling the quasi-linear and non-linear regime 
of structure formation using adapted iV-body codes [11-14], and we are now beginning to 
reach the level of accuracy required for what has been dubbed “precision cosmology”. Such 
accuracy is at the expense of highly intensive computer algorithms which focus on specific 
models and hence on a reduced subset of the full parameter space that needs to be explored. 
We have recently advocated the development of approximate algorithms [15] to explore a 
broader range of parameter space at the expense of accuracy. 

With a modified Wbody code it is then possible to study the observational consequences 
of gravitational screening by generating a wide suite of simulations which in turn can be used 
to design future surveys. Converting such simulations into mock surveys is a crucial step in 
optimising the scientific returns of cosmological experiments such as Euclid [16], LSST [17], 
WFIRST [18] and SKA [19]. 

One particular cosmological probe — the weak lensing of galaxies by intervening gravita¬ 
tional potentials — has been heralded as a particularly powerful test of gravitational physics 
on large scales. It should supply complementary information to probes of the matter density 
field through spectroscopic and photometric galaxy redshift surveys and, in principle, will be 
insensitive to galaxy bias. Current measurements of weak lensing on large scales are not yet 
competitive with other cosmological probes, but the future looks promising with the surveys 
mentioned above. 

In this paper, we describe the algorithms that can be used to produce lensed maps 
in cosmologies with screening. To produce modified gravity simulations we rely on the code 
described in [11] and for lensing we use the GLAMER pipeline [20, 21]. The methods we propose 
here are the starting point for developing a concerted and comprehensive search for signatures 
of screening in observations of weak lensing. Our work complements the work developed in 
[20-25]. 

The outline of this paper is as follows. In Section 2 we briefly summarise the theoretical 
framework with a particular and systematic focus on the different mechanisms for gravita¬ 
tional screening and how they can be parametrised. In Section 3 we describe the modified 
Wbody algorithm with which we produce the density field and the lensing pipeline, MapSim 
and GLAMER. In Section 4 we present the results of a suite of simulations, focusing on maps 
and a few of the main one point and two point statistics. Finally, in Section 5 we discuss 
what we have found as well as the limitations and future applications of our approach. 

2 Theory 

We now establish a useful unified framework for describing the theories of gravity we have 
studied in this paper. 

2.1 Scalar-tensor theories 

In this paper we will focus on scalar-tensor theories of modified gravity that display some 
form of a screening mechanism. These are encompassed by the general Lagrangian^ 

£ = -L d(j), dd(j)) -L V’m), (2.1) 

^ f{R) gravity can be brought to this form by performing a conformal transformation g^i, —> with 

0(p _ 

A{(f>) = e“pi and /3 = see e.g. [26]. 
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where Mpj = (SvrG) ^ and '0m are the matter-fields. Matter is coupled to the scalar field via 
the effective metric 

d/iu — ^ ( 2 - 2 ) 

and therefore moves on geodesics of the metric. In the Einstein frame, where gravity is 
described by general relativity, this implies the existence of a fifth-force; for a test mass in 
the non-relativistic limit the force is given by 

F^ = VlogA(0). (2.3) 

The key to a successful modified gravity model is to have a way of evading the stringent 
constraints from local gravity experiments and at the same time give rise to interesting 
astrophysical and cosmological signatures. One such way is through what is commonly 
referred to as a screening mechanism and below we will give a short review of the different 
screening mechanisms we investigate in this paper. 

2.2 Screening mechanisms 

To see how screening can emerge in a scalar tensor theory let us expand the general scalar 
field Lagrangian Eq. (2.1) to quadratic order about a given field value 0o 

biquadratic = X^"(0o)5;.505.(I0 + (2.4) 

Mpi 

where m{(j)o) = ^(PV{(^q)/ dcj)^ is the (local) value of the scalar field mass, E(0) is the 
potential for 0, and j3 = is the (local) coupling strength of the fifth-force. The 

value 2/3^ = 1 corresponds to a force with the strength of gravity. The fifth-force on a test 
mass, from a source of mass M, can schematically be written in the form 

~ (2.5) 

Let us now consider two different regions of space-time, one where 0o = 0A and another 
where 4>o = <Pb p 4>A respectively. From the quadratic Lagrangian (2.4) we can see that there 
are at least three ways that we can have some form of screening. One way to reduce the effect 
of the fifth-force (in region B compared to A) is by having a large local mass m{4>B) ^ m{4>A) 
which implies a very short interaction range (the chameleon mechanism). If the matter 
coupling l3{(j)B) f3{(j)A), the fifth-force will also be suppressed (the symmetron mechanism). 

Lastly, |A^^(0s)| S> |A^^(0^)| leads to, after canonical normalisation, a weakened matter 
source and therefore also a weakened fifth-force (the Vainshtein mechanism). 

We should note that this simplified description does not tell the whole story; added 
to these effects we can (and do) have additional screening which can only be studied by 
considering the full non-linear dynamics of the scalar field. 

Chameleon mechanism The chameleon mechanism [5, 
by the Lagrangian 

b:0 = ^(v0)2 + v(0). 

The Klein-Gordon equation for the scalar field becomes 

□0 + = 0. (2.7) 


6] can be found in models defined 


( 2 . 6 ) 
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In the presence of matter sources, the dynamics of cj) are determined by an effective potential 
which, for non-relativistic matter, is given by 


Kff = V{(t)) + 


■^(.4^^ Pm 
Mpi 


( 2 . 8 ) 


Too see how screening works in detail, we look at a static, spherically symmetric object of 
density po and radius R embedded in a background of density Penv The approximate solution 
to the Klein-Gordon equation [5] in this case leads to a fifth-force given by 


~ 2/32 


AR GM 


-rrien^r 


(2.9) 


where menv = \/I 4 fr, 00 ( 0 env) is the mass of the scalar field in the background, j3 = 
is the coupling strength (which is constant in the chameleon model) and 


AR 


min 


^ \ 4^0 0env| 1 


( 2 . 10 ) 


is the so-called thin-shell screening factor where d>Ar is the Newtonian potential of the object 
and subscript ”0” (”env”) refer to quantities at the centre (in the environment) of the object. 
The larger the Newtonian potential becomes, the smaller is, and the fifth-force is screened. 
On top of this, and in dense environments, the term \4>o — 4>env\ also becomes small (and ruenv 
becomes large), giving rise to an environmental screening effect, i.e. an object that is not 
screened on its own can be screened if it is in a very dense environment. This environmental 
dependence has many interesting consequences, see e.g. [27-29]. 

Symmetron mechanism The original symmetron model [7] is defined by the same La- 
grangian as for the chameleon where 

^(<(’) =( 2 . 11 ) 

and where p, M, A are model parameters. The effective potential in the presence of matter 
sources is given by 

‘ - 1 ) (2,13) 

The symmetron mechanism is similar to the chameleon mechanism, except that the coupling 
/3 is now field dependent, leading to an additional screening effect. If the local density satisfies 
Pm > the effective potential has a minimum at (/> = 0 near which the field will reside. 

Since the coupling is proportional to cj), the effective matter coupling is suppressed in high 
density regions and the fifth-force is additionally screened. 

The same screening condition Eq. (2.9)-(2.10) as we had for the chameleon also applies 
for the symmetron, but now /3 = /3((^env) is not a constant anymore. 

Vainshtein mechanism The Vainshtein mechanism [10] is responsible for the viability of 
massive gravity, but it can be present in other theories, most notably the Galileons [30] and 
the DGP model [31]. We will not work directly with the original DGP model, but instead 
take a toy model that has the same fifth-force, with the corresponding Vainshtein screening. 
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but where the background is that of ACDM. This is known as the normal branch DGP 
model where we have added dark energy in the form of a cosmological constant [32], The 
field equation for the scalar field, which in the DGP model is the so-called brane-bending 
mode and describes the curvature of the 4D brane we are confined to live in, reads 


□</> + 


3/3DGp(a)a^ 


((□</.)2 - (V^V.</-)2) = 


SttGo? 

3/3dgp 


5 Pn 


(2.14) 


where Tc is the so-called cross-over scale and a is the scale parameter of the background 
metric. In the original DGP model, Tc dictates at what length-scales gravity becomes 5D 
and /3DGp(fl) = 1 + -|- 3 ^^- For static spherical symmetric configurations, the 

field equation reduces to 


1 d 

1 

1 

(M 

2rl d 

1 ^ 

1 

87rG5pm 

r2 

dr 

3/3dgp(«) dr 

\dr ) 

3/3DGp(a) 


This equation can be integrated to yield 

d(j) 2rl ( dct^y _ 2a^ GM 

rdr 3f3-oQp{a)a‘^ \rdr) 3/3DGp(fl) ’ 

which results in the fifth-force on a test-mass being 


GM 1 2 \/l + (rv/r)3 - 1 

3/3DGp(a) (ry/r)3 


(2.15) 


(2.16) 


(2.17) 


where ry = ( so-called Vainshtein radius. The screening mechanism works 

so that the fifth-force is screened at distances smaller than the Vainshtein radius (r <C ry)- 
As opposed to the chameleon and symmetron screening mechanisms, where we have a mass- 
term which gives the fifth-force at finite range, gravity in this case is modified on the largest 
scales: ~ "W sPuUa) 

2.3 Evolution of linear perturbations in modified gravity 

The evolution of linear matter density perturbations in the models we are interested in can 
be written in the general form [33] 

dm + 2 mm = lnm{a)H^6m^^^p^, (2.18) 

where Ges{k,a) is the effective gravitational constant in Fourier space. If Gefr(A:,a) = G we 
recover the equation governing the evolution of the density perturbations in AGDM. 

For the chameleon and symmetron mechanism we have [33] 


Gef[{k,a) _ 2(3‘^{a)P 


(2.19) 


where /3(a) and m{a) are the respective coupling and mass of the scalar field along the 
cosmological attractor solutions. The form of this equation can be understood by noting that 
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it is the ratio of the Fourier transform of the fifth-force potential to the Fourier transform of 
the Newtonian gravitational potential J^[V^<hAr -|- log = 1 -|- 2/3^(a) 

For large scales A:/a <C we have ss l and we recover the ACDM evolution. 

On small scales, k/a ^ l/m{a), = 1 -|- 2/3^ (a) and gravity is modified. For the 

symmetron we have 


m(o) = -^a/I - (ossb/a)^, 
/3(a) = /3o\/l- (ussb/a)^, 


( 2 . 20 ) 

( 2 . 21 ) 


where is the range of the symmetron field at 2 = 0 , /3o = [^^2 i® 

pling strength relative to gravity and Ogsb = is the scale-factor at which the 

modifications of gravity become noticeable. 

For the Hu-Sawicky f{R) model [34], which is one of the models we have simulated, we 


have 


m{a) = 




y/{n+ I)!/™! \ + 4:f2A / 




( 2 . 22 ) 

(2.23) 


where n and |/ro| are model parameters. In this paper we will only consider the case where 
the primordial power spectrum index is n = 1. The DGP model we work with corresponds 
to taking 


/3(a) = 


m{a) = 0 , 

1 


(2.24) 

(2.25) 


\/6/3DGp(a) 

Note that, since m(a) = 0, there is no scale (or k) dependence in Gefr(a). 

Another way to look at these models is within the 7 , /i parameterisation (see e.g. [35]), 
where the metric potentials are given by = AttGo?^ 5pm and 4> = y'k. We have 

m‘^{a)a^ + (1 + 2/3^(a))/c^ 


k- = 


7 = 


m?{a)a? + kk ' 

m?‘{a)a? -|- (1 — 2j3‘^{a))k‘^ 
m?{a)a? -|- (1 -|- 2(3'^{a))k'^' 


(2.26) 

(2.27) 


For the symmetron and f{R) gravity the two functions interpolate between /r = 7 = 1 
for large scales {k/a <C m{a)) to p = 1 + 2/3^(a) and 7 = 1 — small scales 

{k/a ^ m{a)). On the other hand, in our DGP-like model we have p = 1 + 2/3^(a) and 
7 = 1 — i^ 2 / 3 ^°(i) gravity is modified even on the largest scales. 

Linear theory is useful for obtaining a qualitative understanding of what signatures to 
expect, but it neglects an important part of the model’s behaviour, namely the screening 
mechanism. N-body simulations of modified gravity models, such as the ones considered in 
this paper, have shown that the predictions of linear perturbation theory become inaccurate 
as soon as the evolution of the density perturbations starts to deviate from ACDM [36-41]. 
To obtain accurate predictions for the non-linear evolution of these models, we therefore need 
N-body simulations. 
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2.4 Gravitational leasing in modified gravity 

Gravitational leasing is determined by the so-called leasing potential For the 

modified gravity models we consider in this paper this is equivalent to the Newtonian poten¬ 
tial. This can most easily be seen from the 7 , /i parameterisation mentioned in the previous 
section, which gives the following prediction 


$+ = 


(1 + 7)/^ 
2 


^N, 


(2.28) 


where = 4:7rGa‘^5pm is the Newtonian potential, i.e. the same equation as in ACDM. 

When = 2, which is the case for the particular modified gravity theories we study 

in this paper, leasing itself is not modified and the only differences in leasing with respect to 
ACDM are encoded in the differences in the matter distribution caused by the modifications 
of gravity during the process of gravitational collapse. Modified gravity models where leasing 
itself is modified (see e.g. [41]) can be studied within our numerical framework, but we leave 
this for future work. 


3 Simulations 

The numerical results we obtain in this paper have two core elements which we now describe 
in turn: a modified Wbody solver to generate the density field and a weak leasing pipeline 
that can convert a density field into a leasing map. 

3.1 Gravitational Wbody code 

The simulations in this paper have been performed using the ISIS code [11] which is a 
modified gravity modification of the multi-purpose A^-body code RAMSES [42]. The DGP 
model has been implemented following the description laid out in [40]. 

Standard dark matter A^-body simulations are evolved using only two equations. First 
the gravitational potential is calculated (having first used the location of the A^-body particles 
to calculate the density field) using the Poisson equation 

= 47rGa^(pm(3-1) 

and the particles are then evolved using the geodesic equation 

X + 2Hx = (3-2) 

When going to our modified gravity models, the only change (when the background is close 
to AGDM) is that we need to include the fifth-force. This adds a term (see Eq. (2.3)) 
to the right hand side of the geodesic equation. The expression for the force terms and 
the corresponding field equations that we solve in the A^-body code can be found in Ap¬ 
pendix A. Solving these highly non-linear differential equations is the most challenging and 
time-consuming part of modified gravity simulations.^ All the models we have simulated 
have been simulated before; for more information and details about the implementation of 
the scalar field solver and modified gravity simulations in general we refer the reader to 

^See [14] for a comparison of different codes used to simulate modified gravity. The code we used here was 
found to compare very well (to per cent accuracy) with other high-resolution codes deep into the non-linear 
regime. 
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Model 

Type Parameters 

F5 

Fq 

Symm. A 
Symm. B 
DGP, r = 1.2 
DGP, r = 5.6 

~f{R) 

f{R) I/mI = 10-6, n = 1 

Symmetron A<^ = 1.0 Mpc/h, assB = 0.50, /3o = 1-0 

Symmetron = 1.0 Mpc/h, assB = 0.33, /Iq = 1-0 

DGP rcHo = 1.2 

DGP rcHo = 5.6 


Table 1. The modified gravity models for which we have performed fV-body simulations. The 
background cosmology in all the simulations is a standard ACDM cosmology with = 0.271, 
Ha = 0.729, h = 0.703, = 0.966 and as = 0.8. 

[11, 12]. The code we used was also part in a recent code comparison project for modified 
gravity A^-body codes [14] and it was found to agree to the ~ 1% level deep into the non-linear 
regime (A; ~ 5 /i/Mpc) with similar codes for the models simulated for this paper. 

The simulations we have performed all start from the same initial conditions and are 
run with N = 512^ particles in a box of size Bq = 250.0 Mpc//i and with a ACDM cosmology 
given by Clm = 0.271, Da = 0.729, h = 0.703, Ug = 0.966 and cjg = 0.8. For each simulation 
we choose to output the particles needed for the lensing analysis at redshifts z = 0.000, 0.046, 
0.111, 0.176, 0.244, 0.333, 0.422, 0.538, 0.660, 0.818, 1.000, 1.250, 1.500, 1.750, and 1.981. 

The modified gravity models we simulate are the Hu-Sawicky f{R) model, the sym- 
metron model, and the normal branch DGP model. See Table 1 for the parameters used in 
the simulations.^ The background evolution in all the modified gravity simulations is the 
same as in ACDM, allowing for direct comparison of the effect of the fifth-force (and the 
corresponding screening mechanism). 

3.2 Lensing pipeline 

Light passing through an inhomogeneous matter field is deflected by the intervening large- 
scale structure. This effect, called cosmic shear, promises to be a powerful probe of cosmology. 
As long as T = holds for the gravitational potential, the cosmic deflection potential for a 
light cone out to comoving distance x is given by (e.g. [43, 44]) 

where Da is the comoving angular diameter distance. The deflection potential is sourced by 
an effective dimensionless surface mass density 

3H^a„ p Da(x') Da{x - x') s(Da(x') e, x') ^ , 

= - DMx) - 

where 5 is the matter density contrast ^p/p- 

Extracting lensing quantities from simulations using equations (3.3) and (3.4) is cum¬ 
bersome and numerically intensive, as each calculation and each desired source redshift in¬ 
volves raytracing through the simulation volume. Instead, we simulate lensing by large-scale 
structure using the multi-plane approach of the GLAMER lensing pipeline [20, 21]. An ob¬ 
served light cone is first segmented in the radial direction into a number of slices, and the 

^With the parameter choices we have made the DGP model with TcHq = 1.2 (tcHo = 5.6) have the same 
value of ag as the f{B) model with /aq = 10“® {Irq = 10“®). 
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Figure 1. Comparison of a 1.85 deg x 1.25 deg field of view of the convergence k in ACDM and the 
modified gravity models, using a source redshift Zg = 2. All four panels show the same light cone 
extracted from the respective simulations. The convergence maps were simulated with the GLAMER 
pipeline for gravitational leasing. 


three-dimensional mass distribution in each slice is projected onto a plane at the mean co¬ 
moving distance of that slice. The resulting two-dimensional surface mass density maps Ki, 
i = 1,2,... serve as lensing planes for GLAMER, which traces the propagation of light from 
plane to plane using the deflection angle 


OLi = VV’i (3.5) 

and the Poisson equation [21] 

= 2Ki (3.6) 

which relates the surface mass density Ki to the Laplacian of the deflection potential V^V’* 
on each plane. Having thus constructed the lensing simulation, we are free to place a source 
plane at any redshift inside the light cone (i.e. a delta distribution of source redshifts), and 
calculate maps of lensing quantities such as the convergence k or the shear 7 for an observer 
at redshift zero. Sample convergence maps for the ACDM, f{R), symmetron, and DGP 
simulations are shown in figure 1. 

The construction of light cones from simulations and the projection of the mass distri¬ 
bution onto individual planes is done by the MapSim tool [22] in a single step. Each light cone 
is constructed up to redshift Zmax = 2.0, which is thus the highest source redshift available 
for the lensing maps. The field of view of the light cone is a square with a side length of 
3.85 deg, which is the angle subtended by the simulation box size Bq at redshift giving 
us a total area of 14.82 deg^. 

Figure 2 shows the schematic construction of a light cone for the ACDM simulation. It is 
clear that for lower redshifts, much of the simulation box volume is unused. Randomisation 
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z = 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 



comoving distance [Mpc/h] 


Figure 2. Light cone construction with MapSim for the ACDM simulation. The simulation box size is 
Bq = 250 Mpc//i. The dark shaded area is the light cone, growing to a transverse comoving size equal 
to the box size Bq. The mass distribution in each of the individual light cone segments is projected 
onto a leasing plane at the centre of the segment (dashed line). The regions above the light cone 
indicate the simulation snapshots used to construct each segment. The regions below the light cone 
indicate groups of segments that have been randomised in the same way. 


of the box volume through rotation and translation offers a way to extract multiple light 
cones from a single simulation [45], where each of the light cones ends up with a random 
portion of the simulation box. 

MapSim constructs a light cone by randomly picking an origin and orientation of the box 
volume, defining zero comoving distance and the direction of the line of sight. Increasing the 
comoving distance, particles in the field of view are mapped into the light cone, making use 
of the periodic boundary conditions of the simulation box. This continues until the comoving 
distance is a multiple of the simulation box size Bq, changing snapshots as they become a 
better fit for current redshift. The process is repeated, starting from the randomisation, until 
the whole light cone has been constructed. 

Using the technique laid out above, we extract ten randomised light cones from each of 
the simulations and simulate the gravitational lensing of the contained large scale structure. 
We then create lensing maps of the convergence fields k for source redshifts zs = 0.5, 1.0, 
1.5, 2.0. These maps have a size of 2048 x 2048 pixels for the aforementioned field of view of 
3.85 deg x3.85deg, resulting in an angular resolution of 6.77 arcsec per pixel. We note that 
since the simulation box fills the whole field of view at z = 2.0, the same large-scale structure 
is present in each light cone for the higher redshift slices. We thus expect to underestimate 
the sample variance with increasing redshift. 


4 Results 

We now discuss our results, briefly focusing on the power spectrum of the density field before 
we turn to the convergence power spectrum and its evolution. 
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Figure 3. The fractional difference in the matter power spectrum with respect to ACDM for the 
f{R), symmetron, and DGP models. The thin dotted lines shows the linear perturbation theory and 
the dashed lines shows the halofit predictions (using the halofit fitting function of MGCAMB for 
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4.1 Matter power spectrum 

The matter power spectrum is defined via 

{6m{k)6m{q)) = (27r)35(3)(fc + q)P{k), (4.1) 

where 5m{k) is Fourier transform of the density field 6m,{x) = pm{x)/~p^ — 1. We have 
calculated the matter power spectrum for all of our simulations using POWMES [46]. Our 
obtained power spectra agree well with previous Wbody simulations of the same models [ 11 - 
13, 36, 37, 40, 47], which serves as a good cross-check for the validity of our results. We have 
also computed both the linear predictions, using Eq. (2.18), and the halof it [48] predictions 
for the non-linear power spectrum using MGCAMB [49]. This latter code is a modified gravity 
modification of the Boltzmann-code CAMS [50] and uses the improved fitting formula (designed 
for ACDM) of [51] to get better agreement with simulations on small scales. The fractional 
difference of the matter power spectrum with respect to ACDM as a function of redshift can 
be seen in Fig. 3 for the /(R), symmetron, and DGP models, respectively. 

Before we discuss the results in more detail, we note that, overall, the results we find con¬ 
firm the qualitative features we discussed in the theory section. For f{R) and the symmetron 
model, we see that the power spectrum approaches ACDM on large scales, while for DGP 
there are modifications on all scales. Furthermore, comparing the full simulation results with 
the predictions of linear theory shows the effect of screening: linear theory greatly overpre¬ 
dicts the amount of clustering on small scales. This confirms the point we have made above: 
non-linear effects are crucial for accurately modelling the effects of gravitational screening 
and methods, such as the ones presented in this paper, play an essential role. The halof it 
predictions will be discussed in more detail below. 

fiR) For this model, modifications to gravity boost structure formation on small scales. 
Furthermore, we find that P{k) ~ P{k)ACDM for scales k < 0.05 h/Mpc and over all times. 
The effects are stronger in the F 5 model than in the Fq model and this is due to the larger 
range (see Eq. (2.22)) of the fifth-force in the former simulation. For the F 5 {Fq) model we 
see that the matter power spectrum is enhanced by up to Ri 25% (ss 2-5%) at 2 ; = 0 for 
non-linear scales k ~ 1-10 /i/Mpc. At 2 : = 2 the deviations from ACDM are below 5% and 
1% for all scales in the Fq and Fq model respectively. At earlier times the modifications 
are even smaller and thus, for these models, modifications to gravity only have a potentially 
observable impact on structure formation at late times, 2 ; < 2. The main reason for this is 
that the range of the scalar field (again see Eq. (2.22)) decreases rapidly as we go to high 
redshift. Comparing the results of linear perturbation theory with our simulation results, we 
see that the amount of clustering is greatly overpredicted in linear theory. For example, at 
k = 1 /i/Mpc and 2 ; = 0, linear theory predicts a modification of ~ 50% in the Fq model 
whereas the actual result is closer to ss 20 %. 

Symmetron Modifications of gravity in the symmetron model follow the same pattern 
as in f{R) gravity, where structure formation is boosted on small scales. The effects are 
stronger in the B model than in the A model. The reason for this is the smaller value of Ogsb 
in the former simulation, which means that the fifth-force has been active for a longer period 
of cosmic time. The fifth-force is not in operation before the time of symmetry breaking 
= Ossb and consequently the power spectrum is the same as in ACDM for earlier times 
a < Ossb- For the A (B) model we see that the matter power spectrum is enhanced by up to 
~ 15% (~ 5%) at 2 ; = 0 for non-linear scales k ~ 1-10 h/Mpc. Linear theory is an even worse 
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fit to the simulation results in the symmetron model than in f{R) gravity. At fc = 1 /i/Mpc 
at z = 0 we find a ~ 2% modification for the B model whereas linear theory predicts almost 
20%. This larger deviation can be attributed to the stronger screening mechanism in the 
symmetron model, i.e. the fact that oc 4> leads to additional screening on small scales (and 
high density regions) where (p is clustered close to </> = 0. 

DGP Contrary to the two models discussed above, gravity is modified on all scales in the 
DGP model. At z = 0 we find that structure formation in the VcHq = 1.2 (= 5.6) model 
is enhanced by ~ 12% (~ 3%). Going to earlier times the modifications become smaller 
and, at z = 2, are less than « 4% and Ri 1% in the two models, respectively. Since /3dgp 
increases with increasing Tc we expect that a larger Tc leads to stronger modifications, and 
this is indeed the case in our results. On highly non-linear scales k > 5 h/Mpc, deviations 
in the power spectrum are seen to drop for both simulations and at fc = 10 /i/Mpc we are 
very close to the ACDM prediction. This is due to the Vainshtein mechanism being in play 
and reducing the effects of the fifth-force on the small-scale structure formation. We also 
see that for the two models we have chosen to simulate, the relative difference AP/Pacdm 
has a similar shape at all three redshifts depicted in Fig. 3, while the amplitude is markedly 
different. 


4.2 Convergence power spectrum 

Similar to the power spectrum (4.1), the convergence power spectrum Ci is formally defined 
by 

{k{i)k{£')) = (27r)2 + £') Ce , (4.2) 

where £ denotes the angular mode. The convergence power spectrum is an important cos¬ 
mological probe, as it can be related to the evolution history of the universe via Limber’s 
approximation through 


Ciizs) = 


f 


dz 


g‘^{z)il + zf 
X^{z)H{z) 


P{k = llr{z), z) 


(4.3) 


where P{k, z) is the matter power spectrum at redshift z and g{z) is the lensing weight. For 
a single source plane at z = z^ we have g(z} = x(z:)(x(^) ~ x(^s))/x(2s) where x(z) is the 
co-moving distance. This is a powerful combination of many cosmological quantities which 
can help break degeneracies arising in other probes. The convergence power spectrum is also 
directly related to all observed quadratic statistics of cosmic shear measurements, and might 
soon be measured directly (see e.g. [52, 53]). 

In order to estimate the power spectrum from our simulated convergence maps, we use 

C, = J . (4.4) 

where k is the Fourier transform of the convergence field, and the averaging is done over the 
angle of vector £. In practice, we calculate k via a Fast Fourier Transform of the convergence 
maps (after zero-padding to mitigate boundary effects), and perform the averaging in bins of 
A£ with logarithmic spacing. By calculating the convergence power spectrum separately for 
each simulated light cone, we arrive at a sample of results which further gives us a handle 
on the sample variance of our results. The estimated power spectra for our respective f{R), 
symmetron, and DGP models are shown in Fig. 4. 
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Figure 4. Dimensionless convergence power spectra £{£ + 1) for the f{R), symmetron, and DGP 
simulations at source plane redshifts zs = 0.5,1.0,1.5, 2.0 (bottom to top). Shown are the Icr error 
bars over the sample of realisations. Also shown are the predictions obtained from halofit (thick) 
and, in the case of f{R), mghalofit (thin). 









0.001 



Figure 5. We show the relative contribution to the convergence power spectrum from different ranges 
of co-moving Fourier modes of the density fluctuations. The plot is made by computing the halofit 
prediction for Ci and taking P{k, z) = 0 for k values outside the range indicated in the plot. The 
source is here at z = 1.0. 

To estimate up to which i we can trust our results, we computed the relative contribution 
to Ci coming from different ranges of co-moving modes of the matter power spectrum. The 
result can be seen in Fig. 5. The simulations we have performed have a particle Nyquist 
frequency of /^Ny — 6.5 h/Mpc. Fixed-grid simulations with different box sizes usually start 
to deviate from each other for modes larger than k ~ /cNy/4-A;Ny/2, resulting in k^aaiK ~ 
2-3 h/Mpc and £max ~ 2000-3000. However, our simulations have adaptive refinements, 
which means that the effective Nyquist frequency is much larger and a rough estimate from 
the refinement structure gives us a factor 5-10 at 2 ; = 0.0. From this we estimate that we can 
trust the Ci spectra up .^max ~ 10'^. For our largest source redshift Zg = 2.0 this ^max value 
is probably too large due to the lack of refinement at early times. The minimum ^-value we 
can study is fixed by the simulation box size Bq = 250 Mpc//i, corresponding to values of 
fc min = 0.025 /i/Mpc and £min ~ 100 (see Fig. 5). 

4.3 Evolution of the convergence power spectrum 

We will first investigate the redshift evolution of the non-linear power spectrum^ in the 
screened models relative to a reference ACDM simulation (see Fig. 6). The reason for this 
is twofold: First, as noted before, we expect a loss of power in the high z, high i regime 
due to the lack of refinement of the simulations at early times. However, this loss of power 
is consistent among the different simulations, which all evolved with the same refinement 
settings. Although we are deep in the non-linear regime, and the effect of resolution may 

^See [54] for a discussion on the effect of modified gravity on the weak leasing convergence power spectrum 
at linear scales. 
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Figure 6. Left: Dimensionless convergence power spectrum £(£ + 1) for the ACDM simulation 
at source plane redshifts zg = 0.5,1.0,1.5, 2.0 (bottom to top). Shown are the Itr error bars over 
the sample of realisations. Also shown is the prediction obtained from halofit. Right: Difference 
between the ACDM simulation and predictions from CAMB. 
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be different for different models, we expect that the power ratio with respect to 

ACDM will be close to the true value over the full range up to £max ~ and indicative of the 
differences one might find at higher resolution. Another reason for investigating the power 
ratio is of a more physical nature. If is constant in both redshift and mode £, 

it is in principle indistinguishable from ACDM with a different normalisation of the power 
spectrum. If only depends on 2 ; (and not £), it might also be indistinguishable 

from ACDM with a different background history. Thus a robust signature of modified gravity 
requires an evolution in both redshift and scale. 

To robustly estimate Cg /at a given redshift, we calculate the convergence power 
spectra Cg and for the same light cone taken from a screened model and the ACDM 

simulation. Given that all simulations start from the same initial conditions, both power 
spectra trace the evolution of the same patch of the universe, and the ratio Cg/C^^^^ is an 
estimator with much of the influence of cosmic variance taken out. Only after the individual 
ratios Cg/C^^^^ have been calculated in this way do we calculate the mean and variance 
over the sample of light cones. The results for the ratio Cg/C^^^^ can be seen in Fig. 7 for 
the respective f{R), symmetron, and DGP models. As intended, the sample variance has 
been greatly reduced, especially in the region of intermediate £. Below we will discuss the 
results in more details for the three types of models individually. 


fiR) The convergence power spectrum in f{R) shows the same qualitative behaviour as 
the matter power spectrum. On large scales {£ < 0(100)) the fractional difference w.r.t. 
ACDM is close to zero whereas for small angular scales we do find large deviations. In the 
Fq model the deviation from ACDM is consistent (without our error bars) with zero for all 
source redshifts and £ < 10 ^ whereas for the F 5 model we have a non-zero deviation for 
all £ > 10^. The ratio Cg/Cf^^^ increases for both models with the angular scale and at 
I ~ 10^ we find a ^ 20 — 30% (ss 5 — 10%) signal for the F 5 (Fq) model, depending on source 
redshift. 
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Symmetron For the symmetron simulations, we see the same trend in the convergence 
power spectrum as we do for f{R), but with some important differences. The fractional 
difference w.r.t. ACDM grows with i and for I ~ 10“^ we find a « 2-5% (?a 20-30%) deviation 
for the A (B) model depending on source redshift. As opposed to the f{R) model, the 
modifications of gravity are completely absent for a < Ossb for the symmetron. It is only 
close to the present time that the fifth-force has had time to produce a significant difference 
in the matter power spectrum and thus, for the simulation with Ogsb = 0.5, light from high 
redshift sources propagates mostly through a ACDM universe. Because of this, the C^’s for 
the A model are very close to ACDM at z = 2.0. As we move the source redshift closer to the 
present, we see more modifications to large-scale structure, and therefore larger deviations 
in the convergence power spectrum. 

DGP In the DGP model gravity is modified even on the largest scales and we hnd a corre¬ 
sponding deviation at small values of i down to the linear regime. The fractional difference 
w.r.t. ACDM on large scales is found to be of similar order as in the matter power spectrum 
over the range 0 < z < 2. For the VcHq = 1.2 {vcHq = 5.6) model we hnd a ~ 10-15% 
(~ 5-10%) deviation depending on source redshift. The deviation is, within our error bars, 
fairly ^independent, but we do see a slight drop-off of the signal for large that is consistent 
with the results found in the matter power spectrum (see Fig. 3). Contrary to the symmetron 
and f{R) models, we see a redshift-dependency of the signal that is close to a pure scaling 
of the amplitude. 

4.4 Semi-analytical predictions 

We can try and approximate the non-linear effects using semi-analytic methods such as e.g. 
halofit that perform well in the context of ACDM. We hnd that power spectra gener¬ 
ated through a modihed version of halofit can give better agreement with the A^-body 
simulations than linear theory, but even here the predictions for some models (and model 
parameters) can also be very poor. This is to be expected as there is no screening included in 
the formalism and the halofit prediction is a function of the linear P{k, a) only. Neverthe¬ 
less, for specific models, it is possible to obtain good approximations. There do exist other 
approaches for predicting the non-linear matter power spectrum (or equivalently AP/Pacdm) 
than halofit which might do better (see for example [55]), but we have chosen to focus on 
halofit as it is the most widely used method to obtain non-linear predictions for the matter 
power spectrum. Below we will discuss the ability of halofit to predict both the matter 
and convergence power spectra for our three models. 

fiR) For f{R) gravity we saw in Fig. 3 that the linear predictions are a poor match to 
the simulation results for the matter power spectrum and this also holds for the convergence 
power spectrum. The halofit code^ is able to improve significantly on this. For the 
model, the halofit predictions for the convergence power spectrum (see Fig. 4, 7) are roughly 
within the error bars for most of the Grange considered here. Much of this good agreement can 
be attributed to the model-specific fitting formula used in MGCAMB to correct halofit. This 
fitting formula is a 40 parameter fit to f{R) simulations, hence the good level of agreement is 
not surprising. For the Fq model, which is more strongly screened, we see a larger deviation 
for the low redshift sources and for £ = 10^-10^ halofit predicts a ~ 10% signal whereas 
the simulations shows only ~ 5%. 

®halofit in this subsection refers to the modified version provided with MGCAMB [56]. 
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Symmetron For the symmetron there does not exist a fitting formula for the non-linear 
P{k) to be used in halof it. In the bottom panel of Fig. 3 we see that halof it gives a very 
poor fit for the matter power spectrum. It underpredicts the power for scales k < 0.5 /i/Mpc 
and overpredicts the power for k > 0.5 /i/Mpc. However when we use these results to compute 
the convergence power spectrum, see Fig. 4,7, we find a good agreement for the A model 
whereas the halof it predictions show the largest deviations from the simulation results. 
This is a coincidence that comes from the process of integrating the matter power spectrum 
to compute Ci, see Eq. (4.3), for which the over-/underprediction of the power gets averaged 
out. For the B model, we see a larger deviation for the low redshift sources, just as for f{R) 
above, and for i = 10^-10^ halof it predicts a ~ 50% signal whereas the simulations shows 
only Ri 10%. 

DGP For the DGP model, we find that the halof it predictions are fairly good for both 
the matter power spectrum and the convergence power spectrum. The main reason for this is 
that the linear power spectrum is a much better (but far from good enough) approximation 
to the non-linear power spectrum than it is in the other two models considered here. 

In general, the halof it predictions are able to capture the signal, at least qualitatively. 
However, we see a big difference for the model parameters where the screening is largest. 
This illustrates how hard it is to accurately predict the signal of these models on deep non¬ 
linear scales and more elaborate methods are needed to get the accuracy needed to constrain 
these models using data from future weak-lensing surveys in the deep non-linear regime. 


5 Discussion 

In this paper we have investigated the weak leasing signatures of modified gravity theories 
that have a screening mechanism in the deep non-linear regime. We have performed high- 
resolution Wbody simulations of modified gravity models which were selected to cover the 
most common screening mechanisms discussed in the literature. The output of our simula¬ 
tions was then processed by MapSim, which extracts randomised light cones from simulation 
snapshots, collapsing all particles within a given light cone onto a number of planes suitable 
for lensing. These light cones contain the large-scale structure of our modified gravity theo¬ 
ries and, using the GLAMER lensing pipeline, we created realistic convergence maps as observed 
by weak lensing. 

Our analysis of the convergence power spectrum in screened theories of modified gravity 
was performed in two steps: First, we extracted the convergence power spectrum directly 
from the Fourier transform of the generated convergence maps. Our results for ACDM are in 
good agreement with semi-analytical predictions up to £ > 10^, as well as previous studies. 
For higher values of i, we experience a loss in power which we attribute to resolution effects, 
while the lower values of i have a large uncertainty due to the variance within our limited 
sample size. To combat these effects, we have extracted the convergence power spectrum 
relative to a ACDM reference simulation for each individual light cone before 
averaging over realisations. Since all simulations start from the same initial conditions, this 
procedure reduces the influence of randomness by only considering the relative increase in 
power between theories of gravity. The effects of a loss of power due to limited resolution 
are then mitigated and the final result is a clear indicator of the signal we can expect when 
searching for the deviations of screened theories from ACDM. 
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For the particular modified gravity models we have simulated — f{R) gravity, the 
symmetron and the normal branch DGP model — we have chosen parameters that are close 
to the limits that are set by local gravity experiments. For these parameters we have found 
deviations of up to ~ 50% in the convergence power spectrum in the deep non-linear regime 
£ ~ 10^. The scale {£) and source redshift dependence of the signal for these models is 
such that they cannot be mimicked by a ACDM model with a different value of ffm/crg 
which is promising for detecting this signal in future weak leasing surveys. However the deep 
non-linear regime is also where baryonic effects® and the effects of neutrinos on the power 
spectrum are non-negligible. These effects can be highly degenerate with the modified gravity 
signal [58], and a more detailed study is required to quantify these degeneracies. 

We have also compared the simulation results for the matter and convergence power 
spectra to the predictions of linear perturbation theory and those found from applying the 
halof it prescription. Performing modified gravity simulations is computationally expensive, 
so having an accurate prescription for predicting the modified gravity signal in the non-linear 
regime is of great value. We find that linear perturbation theory (meaning the prediction 
for the ratio Pmg / Pacdm is calculated in linear theory) is a very bad approximation to the 
simulation results in almost all cases. This is not surprising since linear theory does not take 
screening into account. The same can be said for halof it, although it is still able to do much 
better than linear theory. Nevertheless, the halof it predictions are still not good enough 
for precision cosmology; one can use detailed fitting formulas such as the one found in the 
MG-halof it code for f{R) gravity, but even there the predictions can be off by as much as 
~ 5-10%. Unfortunately there is no universal fitting formula with this last approach and 
one needs to proceed on a model by model basis; this is clearly unfeasible if one wishes to 
explores extensive swathes of parameter space. 

This paper is a first step in laying down the tools for a thorough analysis of the effects 
of screening in weak leasing. To proceed, we envision a number of steps. For a start, we 
need to develop more efficient methods for generating the realisations of the density field 
(one approach has been advocated in [15]). Indeed multiple realisations will be needed to pin 
down the fine details that will allow us to distinguish between models, and an efficient scan 
of model parameter space is essential to be able to be place reliable constraints on screening 
parameters themselves. 

With such a tool in hand, we need to develop robust analytical methods (in the spirit 
of a modified halof it) which can be incorporated in a likelihood of up and coming data. In 
particular, with the tools we have developed here, it is now possible to calibrate the weak- 
lensing observables of the 3-D simulations with any analytical model we choose to pursue. 
An integral part of this step will be to extend our pipeline to theories where there is non¬ 
trivial gravitational slip (unlike the cases we considered here). This will involve modifying 
the lensing pipeline itself to include the modified integral of the gravitational potentials. 

Finally, it will be possible to focus on specific structure and go beyond the basic statistics 
we have looked at here. As has been shown in [59-61], voids are a promising arena; gravity 
will be unscreened and one expects stronger signatures of modified gravity in such a setting 
(see e.g. [62-66]). 


®See e.g. [57] for a study of baryonic effects on the convergence power spectrum. 
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A Field equations 

In this appendix we present, for completeness, the held equations we solve in the A-body code 
for the different modihed gravity models. In the equations we have applied the quasi-static 
approximation [67-70]. The equation of motion for particles in our A-body simulations can 
be written in the form 


(A.l) 


X + 2Hx = —V<I>Ar — F^, 


where V^<h 7 v = 4:7rGa^{pm — Pm) is the Newtonian potential and is the hfth-force. 
fiR) For the Hu-Sawicky f{R) model we have 



(A.2) 


where the scalar held /r is determined by 



1 j + H^^m 



Symmetron For the symmetron model we have 


Fp = --xVx 



(A.4) 


where the scalar held y is determined by 



(A.5) 


DGP For the DGP model we have 


Fp = 


(A.6) 


where the scalar held (p is determined by 


^2 / j_ ' c 

3/3DGp(a)o^ 


((V2<^)2 _ (V,V,</))2) 


Hi 6pm 

a/^DCP Pm 


(A.7) 


- 21 








References 


[1] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified gravity and cosmology, Physics 
Reports 513 (Mar., 2012) 1-189, [arXiv : 1106.2476], 

[2] J. Khoury, Les Houches Lectures on Physics Beyond the Standard Model of Cosmology, ArXiv 
e-prints (Dec., 2013) [arXiv: 1312.2006]. 

[3] C. M. Will, The Confrontation between General Relativity and Experiment, Living Reviews in 
Relativity 17 (June, 2014) 4, [arXiv: 1403.7377]. 

[4] E. Berti, E. Barausse, V. Cardoso, L. Gualtieri, P. Pani, U. Sperhake, L. C. Stein, N. Wex, 

K. Yagi, T. Baker, C. P. Burgess, F. S. Coelho, D. Doneva, A. De Felice, P. G. Ferreira, 

P. C. C. Freire, J. Healy, C. Herdeiro, M. Horbatsch, B. Kleihaus, A. Klein, K. Kokkotas, 

J. Kunz, P. Laguna, R. N. Lang, T. G. F. Li, T. Littenberg, A. Matas, S. Mirshekari, 

H. Okawa, E. Radu, R. O’Shaughnessy, B. S. Sathyaprakash, C. Van Den Broeck, H. A. 
Winther, H. Witek, M. Emad Aghili, J. Alsing, B. Bolen, L. Bombelli, S. Caudill, L. Chen, 

J. C. Degollado, R. Fujita, C. Gao, D. Gerosa, S. Kamali, H. O. Silva, J. G. Rosa, 

L. Sadeghian, M. Sampaio, H. Sotani, and M. Zilhao, Testing General Relativity with Present 
and Future Astrophysical Observations, ArXiv e-prints (Jan., 2015) [arXiv: 1501.07274]. 

[5] J. Khoury and A. Weltman, Chameleon cosmology, Phys. Rev. D 69 (Feb., 2004) 044026, 
[astro-ph/0309411] . 

[6] D. F. Mota and D. J. Shaw, Evading equivalence principle violations, cosmological, and other 
experimental constraints in scalar field theories with a strong coupling to matter, Phys. Rev. D 
75 (Mar., 2007) 063501, [hep-ph/0608078]. 

[7] K. Hinterbichler and J. Khoury, Screening Long-Range Forces through Local Symmetry 
Restoration, Physical Review Letters 104 (June, 2010) 231301, [arXiv: 1001.4525]. 

[8] M. Pietroni, Dark energy condensation, Phys. Rev. D 72 (Aug., 2005) 043535, 
[astro-ph/0505615]. 

[9] K. A. Olive and M. Pospelov, Environmental dependence of masses and coupling constants, 
Phys. Rev. D 77 (Feb., 2008) 043524, [arXiv:0709.3825]. 

[10] A. 1. Vainshtein, To the problem of nonvanishing gravitation mass. Physics Letters B 39 (May, 
1972) 393-394. 

[11] C. Llinares, D. F. Mota, and H. A. Winther, ISIS: a new N-body cosmological code with scalar 
fields based on RAMSES. Code presentation and application to the shapes of clusters. 
Astronomy and Astrophysics 562 (Feb., 2014) A78, [arXiv : 1307.6748]. 

[12] B. Li, G.-B. Zhao, R. Teyssier, and K. Koyama, ECOSMOG: an Efficient COde for Simulating 
Modified Gravity, JCAPl (Jan., 2012) 51, [arXiv: 1110.1379]. 

[13] E. Puchwein, M. Baldi, and V. Springel, Modified-Gravity-GADGET: a new code for 
cosmological hydrodynamical simulations of modified gravity models, MNRAS 436 (Nov., 2013) 
348-360, [arXiv:1305.2418]. 

[14] H. A. Winther, E. Schmidt, A. Barreira, C. Arnold, S. Bose, G. Llinares, M. Baldi, B. Ealck, 

W. A. Hellwing, K. Koyama, B. Li, D. E. Mota, E. Puchwein, R. Smith, and G.-B. Zhao, 
Modified Gravity N-body Code Comparison Project, ArXiv e-prints (June, 2015) 

[arXiv: 1506.06384]. 

[15] H. A. Winther and P. G. Ferreira, Fast route to nonlinear clustering statistics in modified 
gravity theories, Phys. Rev. D 91 (June, 2015) 123507, [arXiv: 1403.6492]. 

[16] R. Scaramella, Y. Mellier, J. Amiaux, C. Burigana, C. S. Carvalho, J. C. Cuillandre, A. da 
Silva, J. Dinis, A. Derosa, E. Maiorano, P. Pranzetti, B. Garilli, M. Maris, M. Meneghetti, 

I. Tereno, S. Wachter, L. Amendola, M. Cropper, V. Cardone, R. Massey, S. Niemi, 


- 22 


H. Hoekstra, T. Kitching, L. Miller, T. Schrabback, E. Semboloni, A. Taylor, M. Viola, 

T. Maciaszek, A. Ealet, L. Guzzo, K. Jahnke, W. Percival, F. Pasian, M. Sauvage, and the 
Euclid Collaboration, Euclid space mission: a cosmological challenge for the next 15 years, 
ArXiv e-prints (Jan., 2015) [arXiv : 1501. 04908]. 

[17] Z. Ivezic, J. A. Tyson, B. Abel, E. Acosta, R. Allsman, Y. AlSayyad, S. F. Anderson, 

J. Andrew, R. Angel, G. Angeli, R. Ansari, P. Antilogus, K. T. Arndt, P. Astier, E. Aubourg, 
T. Axelrod, D. J. Bard, J. D. Barr, A. Barrau, J. G. Bartlett, B. J. Bauman, S. Beaumont, 

A. C. Becker, J. Becla, C. Beldica, S. Bellavia, G. Blanc, R. D. Blandford, J. S. Bloom, 

J. Bogart, K. Borne, J. F. Bosch, D. Boutigny, W. N. Brandt, M. E. Brown, J. S. Bullock, 

P. Burchat, D. L. Burke, G. Cagnoli, D. Calabrese, S. Chandrasekharan, S. Chesley, E. C. 
Cheu, J. Chiang, C. F. Claver, A. J. Connolly, K. H. Cook, A. Cooray, K. R. Covey, C. Cribbs, 
W. Cui, R. Cutri, G. Daubard, G. Daues, F. Delgado, S. Digel, P. Doherty, R. Dubois, G. P. 
Dubois-Felsmann, J. Durech, M. Eracleous, H. Ferguson, J. Frank, M. Freemon, E. Gangler, 

E. Gawiser, J. C. Geary, P. Gee, M. Geha, R. R. Gibson, D. K. Gilmore, T. Glanzman, 

I. Goodenow, W. J. Gressler, P. Gris, A. Guyonnet, P. A. Hascall, J. Haupt, F. Hernandez, 

C. Hogan, D. Huang, M. E. Huffer, W. R. Innes, S. H. Jacoby, B. Jain, J. Jee, J. G. Jernigan, 

D. Jevremovic, K. Johns, R. L. Jones, C. Juramy-Gilles, M. Juric, S. M. Kahn, J. S. Kalirai, 

N. Kallivayalil, B. Kalmbach, J. P. Kantor, M. M. Kasliwal, R. Kessler, D. Kirkby, L. Knox, 

1. Kotov, V. L. Krabbendam, S. Krughoff, P. Kubanek, J. Kuczewski, S. Kulkarni, R. Lambert, 
L. Le Guillou, D. Levine, M. Liang, K. Lim, C. Lintott, R. H. Lupton, A. Mahabal, 

P. Marshall, S. Marshall, M. May, R. McKercher, M. Migliore, M. Miller, D. J. Mills, D. G. 
Monet, M. Moniez, D. R. Neill, J. Nief, A. Nomerotski, M. Nordby, P. O’Connor, J. Oliver, 

S. S. Olivier, K. Olsen, S. Ortiz, R. E. Owen, R. Pain, J. R. Peterson, C. E. Petry, 

F. Pierfederici, S. Pietrowicz, R. Pike, P. A. Pinto, R. Plante, S. Plate, P. A. Price, M. Prouza, 

V. Radeka, J. Rajagopal, A. Rasmussen, N. Regnault, S. T. Ridgway, S. Ritz, W. Rosing, 

C. Roucelle, M. R. Rumore, S. Russo, A. Saha, B. Sassolas, T. L. Schalk, R. H. Schindler, 

D. P. Schneider, G. Schumacher, J. Sebag, G. H. Sembroski, L. G. Seppala, 1. Shipsey, 

N. Silvestri, J. A. Smith, R. C. Smith, M. A. Strauss, C. W. Stubbs, D. Sweeney, A. Szalay, 

P. Takacs, J. J. Thaler, R. Van Berg, D. Vanden Berk, K. Vetter, F. Virieux, B. Xin, 

L. Walkowicz, C. W. Walter, D. L. Wang, M. Warner, B. Willman, D. Wittman, S. C. Wolff, 

W. M. Wood-Vasey, P. Yoachim, H. Zhan, and for the LSST Collaboration, LSST: from 
Science Drivers to Reference Design and Anticipated Data Products, ArXiv e-prints (May, 

2008) [arXiv: 0805. 2366]. 

[18] D. Spergel, N. Gehrels, J. Breckinridge, M. Donahue, A. Dressier, B. S. Gaudi, T. Greene, 

O. Guyon, C. Hirata, J. Kalirai, N. J. Kasdin, W. Moos, S. Perlmutter, M. Postman, 

B. Rauscher, J. Rhodes, Y. Wang, D. Weinberg, J. Centrella, W. Traub, C. Baltay, J. Colbert, 
D. Bennett, A. Kiessling, B. Macintosh, J. Merten, M. Mortonson, M. Penny, E. Rozo, 

D. Savransky, K. Stapelfeldt, Y. Zu, C. Baker, E. Cheng, D. Content, J. Dooley, M. Foote, 

R. Goullioud, K. Grady, C. Jackson, J. Kruk, M. Levine, M. Melton, C. Peddle, J. Ruffa, and 

S. Shaklan, Wide-Field InfraRed Survey Telescope-Astrophysics Focused Telescope Assets 
WFIRST-AFTA Final Report, ArXiv e-prints (May, 2013) [arXiv: 1305.5422]. 

[19] C. L. Carilli and S. Rawlings, Motivation, key science projects, standards and assumptions, 
NAR 48 (Dec., 2004) 979-984, [astro-ph/0409274]. 

[20] R. B. Metcalf and M. Petkova, GLAMER -LA code for gravitational lensing simulations with 
adaptive mesh refinement, MNRAS 445 (Dec., 2014) 1942-1953, [arXiv: 1312.1128]. 

[21] M. Petkova, R. B. Metcalf, and C. Giocoli, GLAMER - II. Multiple-plane gravitational lensing, 
MNRAS 445 (Dec., 2014) 1954-1966, [arXiv: 1312. 1536]. 

[22] C. Giocoli, R. B. Metcalf, M. Baldi, M. Meneghetti, L. Moscardini, and M. Petkova, 
Disentangling dark sector models using weak lensing statistics, MNRAS 452 (Sept., 2015) 
2757-2772, [arXiv: 1502.03442]. 


- 23 


[23] D. Huterer, Weak lensing and dark energy, Physical Review D 65 (Mar., 2002) 063001, 
[astro-ph/0106399]. 

[24] M. Sato, T. Hamana, R. Takahashi, M. Takada, N. Yoshida, T. Matsubara, and N. Sugiyama, 
Simulations of Wide-Field Weak Lensing Surveys. I. Basic Statistics and Non-Gaussian 
Effects, ApJ 701 (Aug., 2009) 945-954, [arXiv: 0906.2237]. 

[25] A. Cooray, W. Hu, and J. Miralda-Escude, Weak Lensing by Large-Scale Structure: A Dark 
Matter Halo Approach, ApJ Letters 535 (May, 2000) L9-L12, [astro-ph/0003205]. 

[26] P. Brax, C. van de Bruck, A.-C. Davis, and D. J. Shaw, f(R) gravity and chameleon theories, 
Physical Review D 78 (Nov., 2008) 104021, [arXiv:0806.3415]. 

[27] F. Schmidt, Dynamical masses in modified gravity, Physical Review D 81 (May, 2010) 103002, 
[arXiv: 1003. 0409]. 

[28] H. A. Winther, D. F. Mota, and B. Li, Environment Dependence of Dark Matter Halos in 
Symmetron Modified Gravity, ApJ 756 (Sept., 2012) 166, [arXiv: 1110.6438]. 

[29] G.-B. Zhao, B. Li, and K. Koyama, Testing Gravity Using the Environmental Dependence of 
Dark Matter Halos, Physical Review Letters 107 (Aug., 2011) 071303, [arXiv : 1105.0922]. 

[30] A. Nicolis, R. Rattazzi, and E. Trincherini, Galilean as a local modification of gravity, 

Phys. Rev. D 79 (Mar., 2009) 064036, [arXiv:0811.2197]. 

[31] G. Dvali, G. Gabadadze, and M. Porrati, fD gravity on a brane in 5D Minkowski space, 

Physics Letters B 485 (July, 2000) 208-214, [hep-th/0005016]. 

[32] B. Falck, K. Koyama, G.-b. Zhao, and B. Li, The Vainshtein mechanism in the cosmic web, 
JGAP7 (July, 2014) 58, [arXiv: 1404.2206]. 

[33] P. Brax, A.-C. Davis, B. Li, and H. A. Winther, Unified description of screened modified 
gravity. Physical Review D 86 (Aug., 2012) 044015, [arXiv: 1203.4812]. 

[34] W. Hu and 1. Sawicki, Models of f(R) cosmic acceleration that evade solar system tests. 

Physical Review D 76 (Sept., 2007) 064004, [arXiv:0705.1158]. 

[35] T. Baker, P. Ferreira, and C. Skordis, A fast route to modified gravitational growth. Physical 
Review D 89 (Jan., 2014) 024026, [arXiv : 1310.1086]. 

[36] A.-C. Davis, B. Li, D. F. Mota, and H. A. Winther, Structure Pormation in the Symmetron 
Model, ApJ 748 (Mar., 2012) 61, [arXiv: 1108.3081]. 

[37] P. Brax, A.-C. Davis, B. Li, H. A. Winther, and G.-B. Zhao, Systematic simulations of modified 
gravity: symmetron and dilaton models, JGAP 10 (Oct., 2012) 2, [arXiv: 1206.3568]. 

[38] P. Brax, A.-C. Davis, B. Li, H. A. Winther, and G.-B. Zhao, Systematic simulations of modified 
gravity: chameleon models, JGAP 4 (Apr., 2013) 29, [arXiv : 1303.0007]. 

[39] F. Schmidt, Self-consistent cosmological simulations of DGP braneworld gravity. Physical 
Review D 80 (Aug., 2009) 043001, [arXiv:0905.0858]. 

[40] B. Li, G.-B. Zhao, and K. Koyama, Exploring Vainshtein mechanism on adaptively refined 
meshes, JGAP 5 (May, 2013) 23, [arXiv: 1303.0008]. 

[41] A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh, and S. Pascoli, Nonlinear structure 
formation in the cubic Galileon gravity model, JGAP 10 (Oct., 2013) 27, [arXiv : 1306.3219]. 

[42] R. Teyssier, Gosmological hydrodynamics with adaptive mesh refinement. A new high resolution 
code called RAMSES, Astronomy and Astrophysics 385 (Apr., 2002) 337-364, 

[astro-ph/0111367] . 

[43] S. Dodelson, Modern cosmology. Academic Press, 2003. 

[44] P. Schneider, C. S. Kochanek, and J. Wambsganss, Gravitational Lensing: Strong, Weak and 
Micro, vol. 33 of Saas-Fee Advanced Gourses. Springer, Berlin, Heidelberg, 2006. 


- 24 - 


[45] B. Jain, U. Seljak, and S. White, Ray-tracing Simulations of Weak Sensing by Large-Scale 
Structure, ApJ 530 (Feb., 2000) 547-577, [astro-ph/9901191]. 

[46] S. Colombi and D. Novikov, “POWMES: Measuring the Power Spectrum in an N-body 
Simulation.” Astrophysics Source Code Library, Oct., 2011. 

[47] H. Oyaizu, M. Lima, and W. Hu, Nonlinear evolution of f(R) cosmologies. II. Power spectrum, 
Physical Review D 78 (Dec., 2008) 123524, [arXiv:0807.2462]. 

[48] R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. 
Thomas, G. Efstathiou, and H. M. P. Couchman, Stable clustering, the halo model and 
non-linear cosmological power spectra, MNRAS 341 (June, 2003) 1311-1332, 

[astro-ph/0207664] . 

[49] A. Hojjati, L. Pogosian, and G.-B. Zhao, Testing gravity with CAMB and CosmoMC, JCAP 8 
(Aug., 2011) 5, [arXiv: 1106.4543]. 

[50] A. Lewis, A. Challinor, and A. Lasenby, Efficient Computation of Cosmic Microwave 
Background Anisotropies in Closed Friedmann-Robertson-Walker Models, ApJ 538 (Aug., 

2000) 473-476, [astro-ph/9911177]. 

[51] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, Revising the Halofit Model for 
the Nonlinear Matter Power Spectrum, ApJ 761 (Dec., 2012) 152, [arXiv: 1208.2701]. 

[52] J. Harnois-Deraps, S. Vafaei, and L. Van Waerbeke, Cravitational tensing simulations - I. 
Covariance matrices and halo catalogues, MNRAS 426 (Oct., 2012) 1262-1279, 

[arXiv: 1202.2332]. 

[53] M. Kilbinger, Cosmology with cosmic shear observations: a review, Reports on Progress in 
Physics 78 (July, 2015) 086901, [arXiv : 1411.0115]. 

[54] F. Schmidt, Weak leasing probes of modified gravity. Physical Review D 78 (Aug., 2008) 

043002, [arXiv: 0805.4812]. 

[55] P. Brax and P. Valageas, Impact on the power spectrum of screening in modified gravity 
scenarios, Phys. Rev. D 88 (July, 2013) 023527, [arXiv: 1305.5647]. 

[56] G.-B. Zhao, Modeling the Nonlinear Clustering in Modified Gravity Models. I. A Fitting 
Formula for the Matter Power Spectrum of f(R) Gravity, ApJs 211 (Apr., 2014) 23, 

[arXiv: 1312.1291]. 

[57] M. White, Baryons and weak tensing power spectra, Astroparticle Physics 22 (Nov., 2004) 
211-217, [astro-ph/0405593]. 

[58] M. Baldi, F. Villaescusa-Navarro, M. Viel, E. Puchwein, V. Springel, and L. Moscardini, 

Cosmic degeneracies - I. Joint N-body simulations of modified gravity and massive neutrinos, 
MNRAS 440 (May, 2014) 75-88, [arXiv: 1311.2588]. 

[59] P. M. Sutter, G. Lavaux, B. D. Wandelt, D. H. Weinberg, M. S. Warren, and A. Pisani, Voids 
in the SDSS DR9: observations, simulations, and the impact of the survey mask, MNRAS 442 
(Aug., 2014) 3127-3137, [arXiv: 1310.7155]. 

[60] K. Bolejko, C. Glarkson, R. Maartens, D. Bacon, N. Meures, and E. Beynon, Antilensing: The 
Bright Side of Voids, Physical Review Letters 110 (Jan., 2013) 021302, [arXiv : 1209.3142]. 

[61] P. Zivick, P. M. Sutter, B. D. Wandelt, B. Li, and T. Y. Lam, Using cosmic voids to distinguish 
f(R) gravity in future galaxy surveys, MNRAS 451 (Aug., 2015) 4215-4222, 

[arXiv: 1411.5694]. 

[62] Y.-G. Cai, N. Padilla, and B. Li, Testing gravity using cosmic voids, MNRAS 451 (July, 2015) 
1036-1055, [arXiv: 1410.1510]. 

[63] N. Padilla, D. Paz, M. Lares, L. Geccarelli, D. Garcia Lambas, Y.-C. Cai, and B. Li, Void 
Dynamics, ArXiv e-prints (Oct., 2014) [arXiv: 1410.8186]. 


- 25 - 


[64] Y.-C. Cai, N. Padilla, and B. Li, Testing Gravity using Void Profiles, ArXiv e-prints (Oct., 
2014) [arXiv: 1410.8355], 

[65] A. Barreira, M. Cautun, B. Li, C. M. Baugh, and S. Pascoli, Weak lensing by voids in modified 
lensing potentials, J. Cosmology Astropart. Phys. 8 (Aug., 2015) 28, [arXiv: 1505.05809]. 

[66] B. Li, G.-B. Zhao, and K. Koyama, Haloes and voids in f(R) gravity, MNRAS 421 (Apr., 2012) 
3481-3487, [arXiv: 1111. 2602]. 

[67] J. Noller, F. von Braun-Bates, and P. G. Ferreira, Relativistic scalar fields and the quasistatic 
approximation in theories of modified gravity. Physical Review D 89 (Jan., 2014) 023521, 
[arXiv: 1310.3266]. 

[68] S. Bose, W. A. Hellwing, and B. Li, Testing the quasi-static approximation in f(R) gravity 
simulations, JCAP2 (Feb., 2015) 34, [arXiv: 1411.6128]. 

[69] C. Llinares and D. F. Mota, Cosmological simulations of screened modified gravity out of the 
static approximation: Effects on matter distribution. Physical Review D 89 (Apr., 2014) 

084023, [arXiv: 1312.6016]. 

[70] H. A. Winther and P. G. Ferreira, Vainshtein mechanism beyond the quasistatic approximation, 
Phys. Rev. D 92 (Sept., 2015) 064005, [arXiv: 1505.03539]. 


- 26 


